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ABSTRACT 

An approach to assessing the delamination propagation capabilities in commercial 
finite element codes is presented and demonstrated for one code. For this 
investigation, the Double Cantilever Beam (DCB) specimen and the Single Leg 
Bending (SLB) specimen were chosen for full three-dimensional finite element 
simulations. First, benchmark results were created for both specimens. Second, 
starting from an initially straight front, the delamination was allowed to propagate. 
Good agreement between the load-displacement relationship obtained from the 
propagation analysis results and the benchmark results could be achieved by selecting 
the appropriate input parameters. Selecting the appropriate input parameters, 
however, was not straightforward and often required an iterative procedure. 
Qualitatively, the delamination front computed for the DCB specimen did not take 
the shape of a curved front as expected. However, the analysis of the SLB specimen 
yielded a curved front as may be expected from the distribution of the energy release 
rate and the failure index across the width of the specimen. Overall, the results are 
encouraging but further assessment on a structural level is required. 


INTRODUCTION 

One of the most common failure modes for composite structures is delamination 
[1, 2]. The remote loadings applied to composite components are typically resolved 
into interlaminar tension and shear stresses at discontinuities that create mixed-mode 
I, II and III delaminations. To characterize the onset and growth of these 
delaminations, the use of fracture mechanics has become common practice over the 
past two decades [1, 3, 4], The total strain energy release rate, Gt, the mode I 
component due to interlaminar tension, G/, the mode II component due to interlaminar 
sliding shear, G//, and the mode III component, Gm, due to interlaminar scissoring 
shear, as shown in Figure 1, need to be calculated. In order to predict delamination 
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onset or growth for two-dimensional problems, these calculated G components are 
compared to interlaminar fracture toughness properties measured over a range from 
pure mode I loading to pure mode II loading [5-7]. A quasi static mixed-mode fracture 
criterion is determined by plotting the interlaminar fracture toughness, G c , versus the 
mixed-mode ratio, Gi/Gt, determined from data generated using pure mode I Double 
Cantilever Beam (DCB) (Gn/Gf= 0), pure mode II End-Notched Flexure (ENF) 
(G///G]=\ ), and Mixed-Mode Bending (MMB) tests of varying ratios, as shown in 
Figure 2a for T300/914C and Figure 2b for C12K/R6376 [8, 9]. A curve fit of these 
data is performed to determine a mathematical relationship between G, and Gu/Gt. 
[10, 11], Failure is expected when, for a given mixed-mode ratio Gu IGt, the 
calculated total energy release rate, Gt, exceeds the interlaminar fracture toughness, 
G c . An interaction criterion incorporating the scissoring shear (mode III), was recently 
proposed by Reeder [12], The edge-cracked torsion test (ECT) is being considered for 
standardization [13, 14], 

The virtual crack closure technique (VCCT) is widely used for computing energy 
release rates based on results from continuum (2D) and solid (3D) finite element 
analyses and to supply the mode separation required when using the mixed-mode 
fracture criterion [15, 16]. An increased interest in using a fracture mechanics based 
approach to assess the damage tolerance of composite structures in the design phase 
and during certification has also renewed the interest in the virtual crack closure 
technique. The VCCT technique was recently implemented into the commercial finite 
element code ABAQUS* and MD NASTRAN™ [17, 18]. The implementation into 
the commercial finite element code SAMCEF ' [19] is a mix of VCCT and the Virtual 
Crack Extension Method suggested by Parks [20]. As new approaches for analyzing 
composite delamination are incorporated in finite element codes, the need for 
comparison and benchmarking becomes important. 

The objective of this study was to create an approach which allows the assessment 
of delamination propagation capabilities in commercial finite element codes. The 
approach is demonstrated for the commercial finite element code ABAQUS® with 
focus on their implementation of the Virtual Crack Closure Technique (VCCT) [17]. 
For this investigation, the Double Cantilever Beam (DCB) specimen with a 
unidirectional layup and the Single Feg Bending (STB) specimen with multi- 
directional layup (as shown in Figure 3) were chosen for full three-dimensional finite 
element simulations. These specimen configurations were chosen, since a number of 
combined experimental and numerical studies had been performed previously where 
the critical strain energy release rates were evaluated [21, 23-25]. First, benchmark 
results were created using models simulating specimens with different delamination 
lengths. For each delamination length modeled, the load and displacement at the load 
point were monitored. The mixed-mode strain energy release rate components were 
calculated along the delamination front across the width of the specimen. A failure 
index was calculated by correlating the results with the mixed-mode failure criterion 
of the graphite/epoxy material. It was assumed that the delamination propagated when 
the failure index reached unity. Thus, critical loads and critical displacements for 
delamination onset were calculated for each delamination length modeled. These 
critical load/displacement results were used as a benchmark. Second, starting from an 
initially straight front, the delamination was allowed to propagate based on the 
algorithms implemented into VCCT for ABAQUS ®. VCCT control parameters were 
varied to study the effect on the computed load-displacement behavior during 
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propagation. It was assumed that the computed load-displacement relationship should 
closely match the benchmark results established earlier. As a qualitative assessment, 
the shape of the computed delamination fronts were also compared to photographs of 
failed specimens. 


SPECIMEN DESCRIPTION 

For the current numerical investigation, the Double Cantilever Beam (DCB) and 
the Single Leg Bending (SLB) specimens, as shown in Figure 3, were chosen. The 
DCB specimen is used to determine the mode I interlaminar fracture toughness, G/c 
(Gn/Gf= 0) [5], The SLB specimen was introduced for the determination of fracture 
toughness as a function of mixed-mode I/II ratio [21, 22], This test may be performed 
in a standard three-point-bending fixture such as that used for the ENF test. By 
varying the relative thickness of the delaminated regions (ti and ti), a different mixed- 
mode ratio may be achieved. This type of specimen was chosen to study mode 
separation. Previously, a number of combined experimental and numerical studies of 
these specimens had been performed where the critical strain energy release rates 
were evaluated [21, 23-25]. 

In general, mode I, mode II and mixed-mode tests are performed on 
unidirectionally reinforced laminates, which means that delamination growth occurs at 
a [0/0] interface and crack propagation is parallel to the fibers. For the current study, a 
DCB specimen made of T300/1076 graphite/epoxy with a unidirectional layup, [0]24, 
was modeled. Although this unidirectional layup is desired for standard test methods 
to characterize fracture toughness, delamination growth between layers of the same 
orientation will rarely occur in real structures. Previously, combined experimental and 
numerical studies on specimens with multidirectional layups were performed where 
the critical strain energy release rates of various interfaces were evaluated under 
model, mode II and mixed-mode conditions [21, 24], Therefore, a SLB-specimen 
made of C12K/R6376 graphite/epoxy with a multidirectional layup was selected. The 
stacking sequence [±3 0/0/-3 0/0/3 0/0 4 /3 0/0/- 3 0/0/- 3 0/3 0/ 1 - 3 0/3 0/0/3 0/0/ - 

3 0/ 0 4 / 3 0/0/3 0/0/ ±3 0] was designated D±30, where the arrow denotes the location of 
the delamination. The material properties are given in Table I. 


METHODOLOGY 
Fracture Criteria 

Linear elastic fracture mechanics has proven useful for characterizing the onset 
and growth of delaminations in composite laminates [3, 4], A quasi static mixed-mode 
fracture criterion is determined by plotting the interlaminar fracture toughness, G c , 
versus the mixed-mode ratio, Gh/Gt- Typical examples are presented in Figure 2 for 
T300/914C and C12K/R6376 carbon epoxy materials. A fracture criterion was 
suggested by Benzeggah and Kenane [11] using a simple mathematical relationship 
between G c and Gu/Gt 
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( 1 ) 


o. - G i, * (a, k - a k )- . 

In this expression, Gi c and Gu c are the experimentally-determined fracture 
toughness data for mode I and II as shown in Figure 2. The factor r/ was determined 
by a curve fit using the Levenberg-Marquardt algorithm in KaleidaGraph™ graphing 
and data analysis software [26]. Fracture initiation is expected when, for a given 
mixed-mode ratio Gh/Gt, the calculated total energy release rate, Gt, exceeds the 
interlaminar fracture toughness, G c and therefore the failure index Gj/G, is equal or 
greater than unity 


^> 1 . ( 2 ) 

G f 

For three-dimensional analysis, which yields results for the scissoring mode Gui, a 
modified definition is introduced where Gs denotes the sum of the in-plane shearing 
components Gn+Gui [27]. This modification becomes necessary if a mixed- mode 
failure criterion, which accounts for all three modes, is not available. For analyses 
where Gm= 0, this definition is equal to the commonly used definition of the mixed- 
mode ratio, Gn /Gt mentioned above. To determine failure along the delamination 
front, the critical energy release rate G, is calculated using equation (1) with Gu = Gs 
at each point along the delamination front. Subsequently, the failure index Gj/G, is 
determined as above. The modified interaction criterion is an integral part of the 
VCCT for ABAQUS” analysis software and may be selected by the user [17]. 

Recently, Reeder [12] suggested an interaction criterion that is based on the 
fracture criterion suggested by Benzeggah and Kenane but incorporates the mode III 
shear [12] 


G c = G Ic + ( G Ilc - G Ic ) ' 


G II + G III 


+ { G IIIc- G IIc)' 


G II + Gn, 


G „ + G iii 


(3) 


which is also an integral part of the VCCT for ABAQUS R analysis software and may 
be selected by the user [17]. 


Virtual Crack Closure Technique (VCCT) 

A variety of methods are used in the literature to compute the strain energy release 
rate based on results obtained from finite element analysis. For delaminations in 
laminated composite materials where the failure is highly dependent on the mixed- 
mode ratio (as shown in Figure 2), the virtual crack closure technique (VCCT) [15, 
16] has been most widely used for computing energy release rates. VCCT calculations 
using continuum (2D) and solid (3D) finite element analyses provide the mode 
separation required when using the mixed-mode fracture criterion. 

Currently, VCCT for ABAQUS R is an add-on capability to ABAQUS R /Standard 
Versions 6.5, 6.6 and 6.7 that provides a specific implementation of the virtual crack 
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closure technique within ABAQUS®. The plane of delamination in three-dimensional 
analyses is modeled using the existing ABAQUS “/Standard crack propagation 
capability based on the contact pair capability [17]. Beyond simple calculations of the 
mixed-mode strain energy release rates along the delamination front - which was 
studied previously [25] - the implementation also offers a crack propagation capability 
in ABAQUS'*'. Once the energy release rate exceeds the critical strain energy release 
rate (including the user-specified mixed-mode criteria as shown in Figure 2), the node 
at the crack tip is released in the following increment, which allows the crack to 
advance [17]. 

In addition to the mixed-mode fracture criterion, VCCT for ABAQUS R requires, 
additional input for the propagation analysis. If a user specified release tolerance is 
exceeded in an increment ( G-G c )/G c > release tolerance, a cutback operation is 
performed which reduces the time increment. In the new smaller increment, the strain 
energy release rates are recalculated and compared to the user specified cutback 
tolerance. The cutback reduces the degree of overshoot and improves the accuracy of 
the local solution [17]. A release tolerance of 0.2 is suggested in the handbook [17]. 

To help overcome convergence issues during the propagation analysis, 
ABAQUS ® provides: 

• contact stabilization which is applied across only selected contact pairs and 
used to control the motion of two contact pairs while they approach each 
other in multi-body contact. The damping is applied when bonded contact 
pairs debond and move away from each other 

• automatic or static stabilization which is applied to the motion of the entire 
model and is commonly used in models that exhibit statically unstable 
behavior such as buckling 

• viscous regularization which is applied only to nodes on contact pairs that 
have just debonded. The viscous regularization causes the tangent stiffness 
matrix of the softening material to be positive for sufficiently small time 
increments. 

Setting the value of the input parameters correctly is often an iterative procedure, 
which will be discussed later. 


FINITE ELEMENT MODELING 

Typical three-dimensional finite element models of Double Cantilever Beam 
(DCB) and Single Leg Bending (SLB) specimens are shown in Figures 4 and 5. The 
specimens were modeled with solid brick elements C3D8I which had yielded 
excellent results in a previous study [25]. Along the length, all models were divided 
into different sections with different mesh refinement. A refined mesh of length 5 mm 
with 20 elements as used for the DCB specimen is shown in the detail of Figure 4a. 
This section length had been selected in previous studies [23, 25] and was also used 
during the current investigation. Across the width, the model was divided into a center 
section and a refined edge section, i, to capture local edge effects and steep gradients. 
These sections appear as dark areas in the full view of the specimen as shown in 
Figure 4a. The DCB specimen was modeled with six elements through the specimen 
thickness (2/i) as shown in the detail of Figure 4a. This model was used to calculate 
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mode I energy release rates and create the benchmark results discussed later. For 
propagation analyses using VCCT for ABAQUS®, the model with a uniform mesh 
across the width, as shown in Figure 4b, was used to avoid potential problems at the 
transition between the coarse and very fine mesh near the edges of the specimen. 

The plane of delamination was modeled as a discrete discontinuity in the center of 
the specimen. For the analysis with VCCT for ABAQUS®, the models were created as 
separate meshes for the upper and lower part of the specimens with identical nodal 
point coordinates in the plane of delamination [17]. Two surfaces (top and bottom 
surface) were created on the meshes as shown in Figure 4. Additionally, a node set 
was created to identify the intact (bonded nodes) region. 

For the SLB specimen, a model with a uniform mesh across the width was used as 
shown in Figure 5. For modeling convenience, the upper and lower ami were modeled 
similar to the model of the DCB specimen. To model the test correctly only the upper 
ann was supported in the analysis as shown in Figure 5. Two plies on each side of the 
delamination were modeled individually using one element for each ply as shown in 
the detail of Figure 5. Since the delamination occurs at an interface between materials 
with dissimilar properties, care must be exercised in interpreting the values for Gy and 
Gn obtained using the virtual crack closure technique. For interfacial delaminations 
between two differing orthotropic solids, the observed oscillatory singularity at the 
crack tip becomes an issue for small element lengths [28, 29]. Hence, a value of crack 
tip element length, A a, was chosen (approximately three ply thicknesses) in the range 
over which the strain energy release rate components exhibit a reduced sensitivity to 
the value of Act. The adjacent four plies were modeled by one element with material 
properties smeared using the rule of mixtures [30, 31]. This procedure did not 
calculate the full A-B-D stiffness matrix contributions of the plies, however, it 
appeared suitable to enforce a reasonable model size. The adjacent element extended 
over the four 0° plies. The six outermost plies were modeled by one element with 
smeared material properties. 


ANALYSIS 

First, models simulating specimens with different delamination lengths were 
analyzed. For each delamination length modeled, the load and displacement at the 
load point were monitored. The mixed-mode strain energy release rate components 
were calculated along the delamination front across the width of the specimen. A 
failure index was calculated by correlating the results with the mixed-mode failure 
criterion of the graphite/epoxy material. It was assumed that the delamination 
propagated when the failure index reached a value of unity. Thus, critical loads and 
critical displacements for delamination onset were calculated for each delamination 
length modeled. These critical load/displacement results were used as a benchmark. 
Second, starting from an initially straight front, the delamination was allowed to 
propagate based on the algorithm implemented into VCCT for ABAQUS®. Input 
parameters were varied to study the effect on the computed load-displacement 
behavior during propagation. It was assumed that the computed load-displacement 
relationship should closely match the benchmark results established earlier. As a 
qualitative assessment, the shape of the computed delamination fronts were also 
compared to photographs of failed specimens. 
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Creating a Benchmark Solution for DCB specimens 

The computed mode I strain energy release rate values were plotted versus the 
normalized width, v/B, of the specimen as shown in Figure 6a. The results were 
obtained from models shown in Figure 4a for seven different delamination lengths a. 
An opening displacement <5/2=1. Omm was applied to each arm of the model. 
Qualitatively, the mode I strain energy release rate is fairly constant in the center part 
of the specimen and drops progressively towards the edges. This distribution will 
cause the initial straight front to grow into a curved front as explained in detail in the 
literature [32, 33]. As expected, the mode II and mode III strain energy release rates 
were computed to be nearly zero and hence are not shown. Computed mode I strain 
energy release rates decreased with increasing delamination length a. 

The failure index Gt/G c was computed based on a mode I fracture toughness 
G/ c = 170.3 J/m 2 for T300/914C (see Figure 2). The failure index was plotted versus the 
nonnalized width, y/B, of the specimen as shown in Figure 6b. For all delamination 
lengths modeled - except for a = 40 mm - the failure index in the center of the 
specimen (y/B= 0) is above unity (Gt/G c >\). 

For all delamination lengths modeled, the reaction loads P at the location of the 
applied displacement were calculated and plotted versus the applied opening 
displacement <5/2 as shown in Figure 7a. The critical load, P cr n, when the failure index 
in the center of the specimen (y/B= 0) reaches unity ( Gr/G c =\ ), can be calculated based 
on the relationship between load P and the energy release rate G [34], 


P 2 dC P 
2 dA 


(4) 


In equation (4), Cp is the compliance of the specimen and dA is the increase in 
surface area corresponding to an incremental increase in load or displacement at 
fracture. The critical load P cr it and critical displacement <5 a7 ,/2 were calculated for each 
delamination length modeled 


G T 

G~c 





(5) 


and the results were included in the load/displacement plots as shown in Figure 7a. 
The results indicate that, with increasing delamination length, less load is required to 
extend the delamination. This means that the DCB specimen exhibits unstable 
delamination growth under load control. Therefore, prescribed opening displacements 
6/2 were applied in the analysis instead of nodal point loads P to avoid problems with 
numerical stability of the analysis. It was assumed that the critical load/displacement 
results can be used as a benchmark. For the delamination propagation, therefore, the 
load/displacement results obtained from the model of a DCB specimen with an 
initially straight delamination of <7=30 mm length should correspond to the critical 
load/displacement path (in red) in Figure 7a. 



Delamination Propagation in a DCB Specimen using VCCT for ABAQUS® 

The propagation analysis was performed in two steps using the model shown in 
Figure 4b for a delamination length 30 mm. In the first step, a prescribed displacement 
{812= 0.74 mm) was applied in two increments which equaled nearly the critical tip 
opening (<5 cn -/2= 0.75 mm) determined in the analysis above for a delamination length 
of <7=30 mm. In the second step, the total prescribed displacement was increased (<5/2= 
2.8 mm). Automatic incrementation was used with a small increment size at the 
beginning (10‘ 4 of the total increment) and a very small minimum allowed increment 
(10" 18 of the total increment) to reduce the risk of numerical instability and early 
termination of the analysis. The analysis was limited to 1000 increments. Initially, 
analyses were performed without stabilization or viscous regularization. Release 
tolerance values between 0.2 and 0.6 were used. Using the parameters, the analysis 
terminated early prior to advancing the delamination. 

In Figures 7b and 8(a-b), the computed resultant force (load P) at the tip of the 
DCB specimen is plotted versus the applied crack tip opening (d/2) for different input 
parameters which are listed in Table II. For the results shown, the analysis terminated 
when the 1000 increment limit set for the analysis was reached. Several analyses 
terminated early because of convergence problems. To overcome the convergence 
problems, the methods implemented in ABAQUS® were used individually to study 
the effects. For the results plotted in Figure 7b, global stabilization was added to the 
analysis. For a stabilization factor of 2xl0" 5 , the stiffness changed to almost infinity 
once the critical load was reached causing the load to increase sharply (plotted in 
blue). The load increased until a point was reached where the delamination 
propagation started and the load gradually decreased following a zigzag curve with 
local rising and declining segments. The gradual load decrease followed the same 
trend as the benchmark curve (in grey) but is shifted toward higher loads. For a 
stabilization factor of 2x1 0‘ 6 (in green), the same zigzag pattern was observed but the 
average curve was in good agreement with the benchmark result. For a stabilization 
factor of 2xl0" 7 (in red), the average was lower than before but was in good agreement 
with the benchmark result until termination after 550 increments due to convergence 
problems. The results obtained for a stabilization factor of 2x1 0" 8 (in black for a 
release tolerance of 0.2) were on top of the previous result. The rate of convergence 
appeared to be slower since only <5/2= 1.14 mm was applied for 1000 increments 
compared to 6/2= 1.24 mm for a stabilization factor of 2x1 0" 6 and the same release 
tolerance (0.2). Changing the release tolerance also appeared to influence the 
convergence as shown in Table II. For a release tolerance of 0.02, the analysis 
terminated after 1000 increments for 6/2= 1.04 mm. For a release tolerance of 0.002, 
the analysis terminated due to convergence problems after 45 1 increments. Changing 
the release tolerance, however, appeared to have no effect on the overall 
load/displacement behavior or the magnitude of the zigzag pattern. 

For the results plotted in Figure 8a, contact stabilization was added to the analysis. 
For all combinations of stabilization factors and release tolerances, a zigzag pattern 
was observed, where the peak values were in good agreement with the benchmark 
result. The zigzag curve is slightly lower. Decreasing the stabilization factors 
appeared to cause a slower rate of convergence which is either seen by smaller 8/2 for 
the same number of analysis increments or early termination of the analysis as shown 
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in Table II. Changing the release tolerance also appeared to influence the 
convergence. However, it appeared to have no effect on the overall load/displacement 
behavior or the magnitude of the zigzag pattern. 

Viscous regularization was added to the analysis to overcome convergence 
problems. Convergence could not be achieved over a wide range of viscosity 
coefficients when a release tolerance value of 0.2 was used as suggested in reference 
[17]. Subsequently, the release tolerance value was increased. The results where 
convergence was achieved are plotted in Figure 8b. For all combinations of the 
viscosity coefficient and release tolerance, a zigzag pattern was obtained, where the 
peak values were in good agreement with the benchmark result. The average results 
are somewhat lower than the benchmark result. Compared to results obtained from 
analyses with global and contact stabilization, the results obtained with viscous 
regularization appear to have a better rate of convergence since a higher opening 
displacement (<5/2= 1.48 mm) was applied during the analysis for the same number of 
total increments (1000). Decreasing the viscosity coefficient appeared to cause a 
slower rate of convergence which was seen by smaller d/2 values for the same number 
of analysis increments as visible in the plots. Lowering the release tolerance also 
appeared to influence the convergence which was either seen by smaller d/2 for the 
same number of analysis increments as visible in the plots or early termination of the 
analysis as shown in Table II. Changing the release tolerance, however, appeared to 
have no effect on the overall load/displacement behavior or the magnitude of the 
zigzag pattern. 

In summary, good agreement between analysis results and the benchmark could 
be achieved for different release tolerance values in combination with global or 
contact stabilization or viscous regularization. Selecting the appropriate input 
parameters, however, was not straightforward and often required several iterations 
where the parameters had to be changed. All results had a zigzag pattern. 

Besides matching the load displacement behavior of benchmark results, a 
delamination propagation analysis should also yield a delamination front shape that 
is representative of the actual failure. An example of delamination front shapes 
observed by opening a tested DCB specimen are shown in Figure 9a [35]. From the 
initial straight delamination front which is formed by the edge of the Teflon insert, the 
delamination develops into a curved thumb nail shaped front. The front remains 
thumbnail shaped if the test is continued and the delamination continues to grow. 
Delamination propagation computed using the model with a uniform mesh across the 
width (Figure 4b) is shown in Figure 9b at the end of the analysis after 1000 
increments. Plotted on the bottom surface (defined in Figure 4b) are the contours of 
the bond state, where the delaminated section appears in red and the intact (bonded) 
section in blue. The transition between the colors indicates the location of the 
delamination front. The initial straight front was included for clarification. The first 
propagation was observed in the center of the specimen as expected from the 
distribution of the energy release rate (Figure 6a) and the failure index (Figure 6b). 
The front propagated across the width of the specimen until a new straight front was 
reached. Subsequently, the propagation starts again in the center. During the analysis, 
the front never developed into the expected curved thumbnail front, and the analysis 
tenninated with a straight front as shown in Figure 9b. This result is somewhat 
unsatisfactory but may be explained by the fact that the failure index in this particular 
example is nearly constant across about 80% of the width of the specimen as shown in 
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Figure 6b. An even finer mesh may be required to capture the lagging propagation 
near the edge. 


Creating a Benchmark Solution for SLB specimens 

The computed total strain energy release rate values were plotted versus the 
nonnalized width, y/B, of the SLB specimen as shown in Figure 10a. The results were 
obtained from models shown in Figure 5 for twelve different delamination lengths a. 
An arbitrary center deflection w= 2.8 mm was applied as shown in Figure 5. 
Qualitatively, the total energy release rate is fairly constant in the center part of the 
specimen and drops towards the edges. Peaks in the distribution are observed at the 
edges. Computed total strain energy release rates decreased with increasing 
delamination length a. The sum of the shear components G,s = Gn+Gui and the mixed- 
mode ratio Gs /Gt were also calculated for each nodal point along the delamination 
front across the width of the specimen (not shown). 

Using the mixed-mode failure criterion for C12K/R6376 (see Figure 2b), the 
failure index Gt/G c was computed for each node along the delamination front and 
plotted versus the normalized width, y/B, of the specimen as shown in Figure 10b. For 
the center deflection applied, the failure index G-j/G, in the center is well below one. 
The failure index is almost constant in the center of the specimen, drops towards the 
edges and increases again in the immediate vicinity of the edge. To reach Gj/G c =l in 
the center of the specimen (y/B= 0), a critical center deflection, w cnh and 
corresponding critical load P cr it, were calculated using equation (5) for all 
delamination lengths modeled. 

For all delamination lengths modeled, the reaction load P at the location of the 
applied deflection were calculated and plotted versus the applied center deflection, w, 
as shown in Figure 11a. The calculated critical center deflection, w <:nh and 
corresponding critical load values, P C nt, were included in the plots. The results 
indicated that, with increasing delamination length, less load is required to extend the 
delamination. At the same time also, the values of the critical center deflection 
decreased. This means that the SLB specimen exhibits unstable delamination growth 
under load as well as displacement control. From these critical load/displacement 
results, a benchmark solution can be created. To define the benchmark, it is assumed 
that prescribed center deflections are applied in the analysis instead of nodal point 
loads P to minimize problems with numerical stability of the analysis caused by the 
unstable growth. Once the critical center deflection is reached and delamination 
propagation starts, the applied displacement must be held constant over several 
increments while the delamination front is advanced during these increments. Once 
the stable path is reached, the applied center deflection is increased again 
incrementally. For the simulated delamination propagation, therefore, the 
load/displacement results obtained from the model of a SLB specimen with an 
initially straight delamination length of a=34 mm should correspond to the critical 
load/displacement path (in red) as shown in Figure 1 la. 


11 



Delamination Propagation in a SLB Specimen using VCCT for ABAQUS® 

The propagation analysis was performed in two steps using the model shown in 
Figure 5. In the first step, a central deflection (w= 3.1mm) was applied in two 
increments which equaled nearly the critical tip opening (w cn ,= 3.23 mm) determined 
in the analysis above. In the second step, the total prescribed displacement was 
increased (w= 5.0 mm). Automatic time incrementation was used with a small initial 
time increment size ( 1 0‘ 3 ) and a very small minimum allowed time increment (10' 17 ) 
to reduce the risk of numerical instability and early termination of the analysis. The 
analysis was limited to 1000 increments. 

In Figure lib and 12(a-b), the computed resultant force (load P) at the center of 
the SLB specimen is plotted versus the center deflection (vv) for different input 
parameters which are listed in Table II. The analysis terminated before the total 
prescribed center deflection was applied. For the results shown, the analysis 
terminated when the 1000 increment limit set for the analysis was reached. Several 
analyses terminated early because of convergence problems. The results computed 
when global stabilization was used are plotted in Figure 1 lb. For a stabilization factor 
of 2x1 0' 5 , the load increased suddenly at the beginning of the second load step (plotted 
in blue). Then, the load continued to increase on a path with the same stiffness as the 
benchmark but offset to higher loads. The load continued to increase until a point was 
reached where delamination propagation started and the load decreased. The analysis 
was stopped by the user. For a stabilization factor of 2x1 O' 6 (in green), the 
delamination growth started at the critical center deflection. In the beginning, the 
load/displacement path followed the constant deflection branch of the benchmark 
result very well. At the transition between the constant deflection branch and the 
stable propagation branch of the benchmark result, the applied center deflection was 
about 2% higher compared to the benchmark. For the stable path, a zigzag pattern was 
observed but the minimum is in good agreement with the benchmark result. 

The results computed when contact stabilization was used are plotted in 
Figure 12a. For a small stabilization factor (lxl0‘ 6 ) and a release tolerance (0.2) 
suggested in the handbook [17], the load dropped and delamination propagation 
started prior to reaching the critical point of the benchmark solution (plotted in blue). 
The load/displacement path then ran parallel to the constant deflection branch of the 
benchmark result but the analysis terminated early due to convergence problems. The 
stabilization factor and release tolerance had to be increased to avoid premature 
termination of the analysis. For a stabilization factor of lxl 0' 5 and release tolerance of 
0.5 (in green), the load dropped at the critical point of the benchmark solution. First, 
the center deflection kept increasing with decreasing load. Later, the 
load/displacement path ran parallel to the constant deflection branch of the benchmark 
result. At the transition between the constant deflection branch and the stable 
propagation branch of the benchmark result, the applied center deflection was about 
2% higher compared to the benchmark. For the stable path, a zigzag pattern was 
observed where the average results were in good agreement with the benchmark 
result. The difference between the maximum and minimum values was much smaller 
than in the case where global stabilization was used. The best results compared to the 
benchmark were obtained for even higher values of the stabilization factors of lxl 0' 4 
and a release tolerance of 0.5 (in red). 
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When viscous regularization was used to help overcome convergence issues, a 
value of 0.2 was used initially for the release tolerance as suggested in the handbook 
[17]. Convergence could not be achieved which led to an increase in release tolerance. 
The results are plotted in Figure 12b. For a small viscosity coefficient of lxl 0" 4 and a 
release tolerance of 0.5 (in blue), the load dropped at the critical point, but the center 
deflection kept increasing with decreasing load. Then, the analysis terminated early 
due to convergence problems. For an increased viscosity coefficient of 1x10 " and a 
release tolerance of 0.5 (in red), the load dropped at the critical point and the 
load/displacement path started following the constant deflection branch of the 
benchmark result, but the analysis terminated early due to convergence problems. The 
viscosity coefficient and release tolerance had to be increased further to avoid 
premature termination of the analysis. For a viscosity coefficient of lxl 0" 1 and a 
release tolerance of 0.9 (in green), the load dropped at the critical point. First, the 
center deflection kept increasing with decreasing load. Later, the load/displacement 
path ran parallel to the constant deflection branch of the benchmark result. At the 
transition between the constant deflection branch and the stable propagation branch of 
the benchmark result, the applied center deflection is about 2.5% higher compared to 
the benchmark. For the stable path, a zigzag pattern is observed where the average 
results are in good agreement with the benchmark result. The difference between the 
maximum and minimum values is much smaller compared to the cases where global 
or contact stabilization was used. 

In summary, good agreement between analysis results and the benchmark could 
be achieved for different release tolerance values in combination with global or 
contact stabilization or viscous regularization. Selecting the appropriate input 
parameters, however, was not straightforward and often required several iterations 
where the parameters had to be changed. 

Delamination propagation computed using the model with a uniform mesh across 
the width (Figure 5) is shown in Figure 13 after 1000 increments. Plotted on the 
bottom surface (defined in Figure 5) are the contours of the bond state variable. The 
bond state varies between 0.0 (fully bonded shown in dark blue) and 1.0 (fully 
disbonded shown in red) [17]. The transition between the colors indicated the location 
of the delamination front. The initial straight front was included for clarification. The 
first propagation is observed near the center and corresponds to the maximum in the 
distribution of the failure index (Figure 10b). The front then propagated across the 
width. Further propagation created a curved front where the edges lag behind as 
shown in Figure 13. This result is in good agreement with expectations based on the 
distribution of the failure index shown in Figure 10b. C-scans or x-ray photographs of 
tested specimens were not available for comparison. 


SUMMARY AND CONCLUSIONS 

An approach for assessing the delamination propagation capabilities in 
commercial finite element codes is presented and demonstrated for the commercial 
finite element code ABAQUS R with focus on their implementation of the Virtual 
Crack Closure Technique (VCCT). For this investigation, the Double Cantilever 
Beam (DCB) specimen with a unidirectional layup and the Single Leg Bending (SLB) 
specimen with a multi-directional layup were chosen for full three-dimensional finite 
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element simulations. First, critical load/displacement results were defined for 
delamination onset which were used subsequently as benchmarks. Second, starting 
from an initially straight front, the delamination was allowed to propagate based on 
the algorithms implemented into VCCT for ABAQUS®. VCCT control parameters 
were varied to study the effect on the computed load-displacement behavior during 
propagation. It was assumed that for good results the computed load-displacement 
relationship should correspond to the benchmark results established earlier. Third - as 
a qualitative assessment - the shape of the computed delamination fronts were also 
compared to photographs of failed specimens. 

Good agreement between the load-displacement relationship obtained from the 
propagation analysis results and the benchmark results could be achieved by selecting 
the appropriate input parameters, however, a zigzag response was obtained during 
propagation. Selecting the appropriate VCCT input parameters such as release 
tolerance, global or contact stabilization and viscous regularization, however, was 
not straightforward and often required an iterative procedure. In this case, the input 
parameters were modified until the analysis results agreed with the benchmark. For 
all the combinations of input parameters, only a global stabilization factor of 2xl0‘ 6 
in combination with a release tolerance of 0.2 gave good results for the DCB and SLB 
simulations. In a real case scenario where the results are unknown, obtaining the 
right solution will remain challenging. 

Besides matching the load displacement behavior of the benchmark results, a 
delamination propagation analysis should also yield a delamination front shape that is 
representative of the actual failure. During the analysis of the DCB specimen, the 
front never developed into the expected curved thumbnail front as seen in tested 
specimens. The analysis terminated with a straight front which is somewhat 
unsatisfactory. The result may be explained by the fact that the failure index in this 
particular example is constant across about 80% of the width of the specimen and a 
finer mesh may be required to capture the lagging propagation near the edge. During 
the analysis of the SLB specimen, the front developed into a curved front as expected 
from the distribution of the failure index. This result is encouraging. Overall, the 
results are promising but further studies are required which should include different 
levels of mesh refinement, new stabilization options and the use of continuum shell 
elements to model the specimens. Additionally, assessment of the propagation 
capabilities in more complex specimens and on a structural level is required. 
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TABLE 1. MATERIAL PROPERTIES. 


T300/1076 Unidirectional Graphite/Epoxy Prepreg 

E\ i = 139.4 GPa 

e 22 = 10.16 GPa 

£33 = 10.16 GPa 

vi 2 = 0.30 

V 13 = 0.30 

V23 = 0.436 

G 12 = 4.6 GPa 

Gl 3 = 4.6 GPa 

G 23 =3.54 GPa 

C12K/R6376 Unidirectional Graphite/Epoxy Prepreg 

E\ i = 146.9 GPa 

£22 = 10 6 GPa 

£33 = 10.6 GPa 

vi 2 = 0.33 

vi 3 = 0.33 

V 23 = 0.33 

Gl 2 = 5.45 GPa 

Gl 3 = 5.45 GPa 

G 23 = 3.99 GPa 


The material properties are given with reference to the ply coordinate axes where index 1 1 denotes 
the ply principal axis that coincides with the direction of maximum in-plane Young’s modulus (fiber 
direction). Index 22 denotes the direction transverse to the fiber in the plane of the lamina and index 
33 the direction perpendicular to the plane of the lamina. 


TABLE II. INPUT PARAMETERS. 


FE model 

global 

stabilization 

contact 

stabilization 

viscous 

regularization 

release 

tolerance 

last increment 

DCB-st3 

2 10 -3 



0.2 

381 

DCB-st4 

2 10 -6 



0.2 

1002 

DCB-st5 

2 10' v 



0.2 

550 

DCB-st6 

2 10 -8 



0.2 

1002 

DCB-st7 

2 10 -8 



0.02 

1002 

DCB-st8 

2 10 -8 



0.002 

451 

DCB-ctl 


Hit 5 


0.2 

1002 

DCB-ct2 


1 10 -6 


0.2 

1002 

DCB-ct3 


1 10 -7 


0.2 

751 

DCB-ct4 


1 10 -7 


0.02 

1002 

DCB-ct5 


1 10 -7 


0.002 

1002 

DCB-ct6 


Hit 3 


0.002 

911 

DCB-vrl 



no^ 

0.5 

1002 

DCB-vr2 



no^ 

0.3 

273 

DCB-vr3 



1 10' 5 

0.5 

1002 

DCB-vr4 



no 15 

0.3 

1002 

SLB-stl 

2 10 -3 



0.2 

266 

SLB-st2 

2 10 -6 



0.2 

1002 

SLB-ctl 


1 10 -6 


0.2 

133 

SLB-ct6 


Hit 5 


0.5 

811 

SLB-ct8 


Ho^ 


0.5 

1002 

SLB-vrl 



no^ 

0.5 

65 

SLB-vr6 



no 12 

0.5 

88 

SLB-vrl 2 



no 71 

0.9 

537 
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Figure 1: Fracture Modes. 
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Figure 2a. Mixed-mode fracture criterion for T300/914C. 
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Figure 2b. Mixed-mode fracture criterion for C12K/R6376. 
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(a) Double Cantilever Beam Specimen (DCB) 
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D30: Cl 2K/R6376 [±30/0/-30/0/30/0 4 /30/0/-30/0/-30/30//-30/30/0/30/0/ -30/0 4 /30/0/30/0/±30] 

(b) Single Leg Bending Specimen (SLB) 

Figure 3. Specimen configurations. 
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a. Deformed model ofDCB specimen with refined edges 



b. Deformed model of a DCB specimen for VCCT for ABAQUS analysis 
Figure 4. Full three-dimensional finite element model of a DCB specimen. 
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Figure 5. Deformed model ofSLB specimen and detail of region around delamination front. 
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Figure 6a. Computed strain energy release rate distribution across the width of 
a DCB specimen (model Figure 4a). 
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Figure 6b. Failure index distribution across the width of a DCB specimen) model Figure 4a). 
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Figure 7a. Benchmark: Critical load-clisplacement behavior for DCB specimen. 



applied opening displacement 8/2, mm 

Figure 7b. VCCT for ABAQUS: Computed critical load-displacement behavior 
for DCB specimen obtained from results with global stabilization. 
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applied opening displacement 5/2, mm 

Figure 8a. Computed critical load-displacement behavior for DCB specimen 
obtained from results with contact stabilization. 



applied opening displacement 6/2, mm 

Figure 8b. Computed critical load-displacement behavior for DCB specimen 
obtained from results with viscous regularization. 
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(a) Scan of fractured DCB specimen 



(b) Delamination front after 1000 increments (Bottom surface ofFE model in Figure 4b) 
Figure 9. Delamination front shape for a DCB specimen. 
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Figure 10a. Computed strain energy release rate distribution 
across the width of a SLB specimen. 
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Figure 10b. Failure index distribution across the width of a SLB specimen. 
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Figure 11a. Critical load-displacement behavior for SLB specimen. 



applied center deflection w, mm 

Figure lib. Computed critical load-displacement behavior for SLB specimen 
obtained from results with global stabilization. 



Figure 12a. Computed critical load-displacement behavior for SLB specimen 
obtained from results with contact stabilization. 
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Figure 12b. Computed critical load-displacement behavior for SLB specimen 
obtained from results with viscous regularization. 



location of the initial delamination front 


Figure 13. Delamination front for a SLB specimen (Bottom surface ofFE model in Figure 5). 


